/*
 * Copyright 2013 Google Inc.
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *      http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */

package com.aplink.generic.google.maps;

import java.util.ArrayList;
import java.util.List;

import com.google.android.gms.maps.model.LatLng;

public class PolyUtil {

    private static final double DEFAULT_TOLERANCE = 0.1; // meters.

    /**
     * Computes whether the given point lies inside the specified polygon. The
     * polygon is always cosidered closed, regardless of whether the last point
     * equals the first or not. Inside is defined as not containing the South
     * Pole -- the South Pole is always outside. The polygon is formed of great
     * circle segments if geodesic is true, and of rhumb (loxodromic) segments
     * otherwise.
     */
    public static boolean containsLocation(final LatLng point,
            final List<LatLng> polygon, final boolean geodesic) {
        final int size = polygon.size();
        if (size == 0) {
            return false;
        }
        final double lat3 = Math.toRadians(point.latitude);
        final double lng3 = Math.toRadians(point.longitude);
        final LatLng prev = polygon.get(size - 1);
        double lat1 = Math.toRadians(prev.latitude);
        double lng1 = Math.toRadians(prev.longitude);
        int nIntersect = 0;
        for (final LatLng point2 : polygon) {
            final double dLng3 = MathUtil.wrap(lng3 - lng1, -Math.PI, Math.PI);
            // Special case: point equal to vertex is inside.
            if ((lat3 == lat1) && (dLng3 == 0)) {
                return true;
            }
            final double lat2 = Math.toRadians(point2.latitude);
            final double lng2 = Math.toRadians(point2.longitude);
            // Offset longitudes by -lng1.
            if (PolyUtil.intersects(lat1, lat2,
                    MathUtil.wrap(lng2 - lng1, -Math.PI, Math.PI), lat3, dLng3,
                    geodesic)) {
                ++nIntersect;
            }
            lat1 = lat2;
            lng1 = lng2;
        }
        return (nIntersect & 1) != 0;
    }

    /**
     * Decodes an encoded path string into a sequence of LatLngs.
     */
    public static List<LatLng> decode(final String encodedPath) {
        final int len = encodedPath.length();

        // For speed we preallocate to an upper bound on the final length, then
        // truncate the array before returning.
        final List<LatLng> path = new ArrayList<LatLng>();
        int index = 0;
        int lat = 0;
        int lng = 0;

        while (index < len) {
            int result = 1;
            int shift = 0;
            int b;
            do {
                b = encodedPath.charAt(index++) - 63 - 1;
                result += b << shift;
                shift += 5;
            } while (b >= 0x1f);
            lat += (result & 1) != 0 ? ~(result >> 1) : (result >> 1);

            result = 1;
            shift = 0;
            do {
                b = encodedPath.charAt(index++) - 63 - 1;
                result += b << shift;
                shift += 5;
            } while (b >= 0x1f);
            lng += (result & 1) != 0 ? ~(result >> 1) : (result >> 1);

            path.add(new LatLng(lat * 1e-5, lng * 1e-5));
        }

        return path;
    }

    /**
     * Encodes a sequence of LatLngs into an encoded path string.
     */
    public static String encode(final List<LatLng> path) {
        long lastLat = 0;
        long lastLng = 0;

        final StringBuffer result = new StringBuffer();

        for (final LatLng point : path) {
            final long lat = Math.round(point.latitude * 1e5);
            final long lng = Math.round(point.longitude * 1e5);

            final long dLat = lat - lastLat;
            final long dLng = lng - lastLng;

            PolyUtil.encode(dLat, result);
            PolyUtil.encode(dLng, result);

            lastLat = lat;
            lastLng = lng;
        }
        return result.toString();
    }

    private static void encode(long v, final StringBuffer result) {
        v = v < 0 ? ~(v << 1) : v << 1;
        while (v >= 0x20) {
            result.append(Character.toChars((int) ((0x20 | (v & 0x1f)) + 63)));
            v >>= 5;
        }
        result.append(Character.toChars((int) (v + 63)));
    }

    /**
     * Computes whether the vertical segment (lat3, lng3) to South Pole
     * intersects the segment (lat1, lng1) to (lat2, lng2). Longitudes are
     * offset by -lng1; the implicit lng1 becomes 0.
     */
    private static boolean intersects(final double lat1, final double lat2,
            final double lng2, final double lat3, final double lng3,
            final boolean geodesic) {
        // Both ends on the same side of lng3.
        if (((lng3 >= 0) && (lng3 >= lng2)) || ((lng3 < 0) && (lng3 < lng2))) {
            return false;
        }
        // Point is South Pole.
        if (lat3 <= (-Math.PI / 2)) {
            return false;
        }
        // Any segment end is a pole.
        if ((lat1 <= (-Math.PI / 2)) || (lat2 <= (-Math.PI / 2))
                || (lat1 >= (Math.PI / 2)) || (lat2 >= (Math.PI / 2))) {
            return false;
        }
        if (lng2 <= -Math.PI) {
            return false;
        }
        final double linearLat = ((lat1 * (lng2 - lng3)) + (lat2 * lng3))
                / lng2;
        // Northern hemisphere and point under lat-lng line.
        if ((lat1 >= 0) && (lat2 >= 0) && (lat3 < linearLat)) {
            return false;
        }
        // Southern hemisphere and point above lat-lng line.
        if ((lat1 <= 0) && (lat2 <= 0) && (lat3 >= linearLat)) {
            return true;
        }
        // North Pole.
        if (lat3 >= (Math.PI / 2)) {
            return true;
        }
        // Compare lat3 with latitude on the GC/Rhumb segment corresponding to
        // lng3.
        // Compare through a strictly-increasing function (tan() or mercator())
        // as convenient.
        return geodesic ? Math.tan(lat3) >= PolyUtil.tanLatGC(lat1, lat2, lng2,
                lng3) : MathUtil.mercator(lat3) >= PolyUtil.mercatorLatRhumb(
                lat1, lat2, lng2, lng3);
    }

    /**
     * Same as {@link #isLocationOnEdge(LatLng, List, boolean, double)} with a
     * default tolerance of 0.1 meters.
     */
    public static boolean isLocationOnEdge(final LatLng point,
            final List<LatLng> polygon, final boolean geodesic) {
        return PolyUtil.isLocationOnEdge(point, polygon, geodesic,
                PolyUtil.DEFAULT_TOLERANCE);
    }

    /**
     * Computes whether the given point lies on or near the edge of a polygon,
     * within a specified tolerance in meters. The polygon edge is composed of
     * great circle segments if geodesic is true, and of Rhumb segments
     * otherwise. The polygon edge is implicitly closed -- the closing segment
     * between the first point and the last point is included.
     */
    public static boolean isLocationOnEdge(final LatLng point,
            final List<LatLng> polygon, final boolean geodesic,
            final double tolerance) {
        return PolyUtil.isLocationOnEdgeOrPath(point, polygon, true, geodesic,
                tolerance);
    }

    private static boolean isLocationOnEdgeOrPath(final LatLng point,
            final List<LatLng> poly, final boolean closed,
            final boolean geodesic, final double toleranceEarth) {
        final int size = poly.size();
        if (size == 0) {
            return false;
        }
        final double tolerance = toleranceEarth / MathUtil.EARTH_RADIUS;
        final double havTolerance = MathUtil.hav(tolerance);
        final double lat3 = Math.toRadians(point.latitude);
        final double lng3 = Math.toRadians(point.longitude);
        final LatLng prev = poly.get(closed ? size - 1 : 0);
        double lat1 = Math.toRadians(prev.latitude);
        double lng1 = Math.toRadians(prev.longitude);
        if (geodesic) {
            for (final LatLng point2 : poly) {
                final double lat2 = Math.toRadians(point2.latitude);
                final double lng2 = Math.toRadians(point2.longitude);
                if (PolyUtil.isOnSegmentGC(lat1, lng1, lat2, lng2, lat3, lng3,
                        havTolerance)) {
                    return true;
                }
                lat1 = lat2;
                lng1 = lng2;
            }
        } else {
            // We project the points to mercator space, where the Rhumb segment
            // is a straight line,
            // and compute the geodesic distance between point3 and the closest
            // point on the
            // segment. This method is an approximation, because it uses
            // "closest" in mercator
            // space which is not "closest" on the sphere -- but the error is
            // small because
            // "tolerance" is small.
            final double minAcceptable = lat3 - tolerance;
            final double maxAcceptable = lat3 + tolerance;
            double y1 = MathUtil.mercator(lat1);
            final double y3 = MathUtil.mercator(lat3);
            final double[] xTry = new double[3];
            for (final LatLng point2 : poly) {
                final double lat2 = Math.toRadians(point2.latitude);
                final double y2 = MathUtil.mercator(lat2);
                final double lng2 = Math.toRadians(point2.longitude);
                if ((Math.max(lat1, lat2) >= minAcceptable)
                        && (Math.min(lat1, lat2) <= maxAcceptable)) {
                    // We offset longitudes by -lng1; the implicit x1 is 0.
                    final double x2 = MathUtil.wrap(lng2 - lng1, -Math.PI,
                            Math.PI);
                    final double x3Base = MathUtil.wrap(lng3 - lng1, -Math.PI,
                            Math.PI);
                    xTry[0] = x3Base;
                    // Also explore wrapping of x3Base around the world in both
                    // directions.
                    xTry[1] = x3Base + (2 * Math.PI);
                    xTry[2] = x3Base - (2 * Math.PI);
                    for (final double x3 : xTry) {
                        final double dy = y2 - y1;
                        final double len2 = (x2 * x2) + (dy * dy);
                        final double t = len2 <= 0 ? 0 : MathUtil.clamp(
                                ((x3 * x2) + ((y3 - y1) * dy)) / len2, 0, 1);
                        final double xClosest = t * x2;
                        final double yClosest = y1 + (t * dy);
                        final double latClosest = MathUtil
                                .inverseMercator(yClosest);
                        final double havDist = MathUtil.havDistance(lat3,
                                latClosest, x3 - xClosest);
                        if (havDist < havTolerance) {
                            return true;
                        }
                    }
                }
                lat1 = lat2;
                lng1 = lng2;
                y1 = y2;
            }
        }
        return false;
    }

    /**
     * Same as {@link #isLocationOnPath(LatLng, List, boolean, double)}
     * 
     * with a default tolerance of 0.1 meters.
     */
    public static boolean isLocationOnPath(final LatLng point,
            final List<LatLng> polyline, final boolean geodesic) {
        return PolyUtil.isLocationOnPath(point, polyline, geodesic,
                PolyUtil.DEFAULT_TOLERANCE);
    }

    /**
     * Computes whether the given point lies on or near a polyline, within a
     * specified tolerance in meters. The polyline is composed of great circle
     * segments if geodesic is true, and of Rhumb segments otherwise. The
     * polyline is not closed -- the closing segment between the first point and
     * the last point is not included.
     */
    public static boolean isLocationOnPath(final LatLng point,
            final List<LatLng> polyline, final boolean geodesic,
            final double tolerance) {
        return PolyUtil.isLocationOnEdgeOrPath(point, polyline, false,
                geodesic, tolerance);
    }

    private static boolean isOnSegmentGC(final double lat1, final double lng1,
            final double lat2, final double lng2, final double lat3,
            final double lng3, final double havTolerance) {
        final double havDist13 = MathUtil.havDistance(lat1, lat3, lng1 - lng3);
        if (havDist13 <= havTolerance) {
            return true;
        }
        final double havDist23 = MathUtil.havDistance(lat2, lat3, lng2 - lng3);
        if (havDist23 <= havTolerance) {
            return true;
        }
        final double sinBearing = PolyUtil.sinDeltaBearing(lat1, lng1, lat2,
                lng2, lat3, lng3);
        final double sinDist13 = MathUtil.sinFromHav(havDist13);
        final double havCrossTrack = MathUtil
                .havFromSin(sinDist13 * sinBearing);
        if (havCrossTrack > havTolerance) {
            return false;
        }
        final double havDist12 = MathUtil.havDistance(lat1, lat2, lng1 - lng2);
        final double term = havDist12 + (havCrossTrack * (1 - (2 * havDist12)));
        if ((havDist13 > term) || (havDist23 > term)) {
            return false;
        }
        if (havDist12 < 0.74) {
            return true;
        }
        final double cosCrossTrack = 1 - (2 * havCrossTrack);
        final double havAlongTrack13 = (havDist13 - havCrossTrack)
                / cosCrossTrack;
        final double havAlongTrack23 = (havDist23 - havCrossTrack)
                / cosCrossTrack;
        final double sinSumAlongTrack = MathUtil.sinSumFromHav(havAlongTrack13,
                havAlongTrack23);
        return sinSumAlongTrack > 0; // Compare with half-circle == PI using
                                     // sign of sin().
    }

    /**
     * Returns mercator(latitude-at-lng3) on the Rhumb line (lat1, lng1) to
     * (lat2, lng2). lng1==0.
     */
    private static double mercatorLatRhumb(final double lat1,
            final double lat2, final double lng2, final double lng3) {
        return ((MathUtil.mercator(lat1) * (lng2 - lng3)) + (MathUtil
                .mercator(lat2) * lng3)) / lng2;
    }

    /**
     * Returns sin(initial bearing from (lat1,lng1) to (lat3,lng3) minus initial
     * bearing from (lat1, lng1) to (lat2,lng2)).
     */
    private static double sinDeltaBearing(final double lat1, final double lng1,
            final double lat2, final double lng2, final double lat3,
            final double lng3) {
        final double sinLat1 = Math.sin(lat1);
        final double cosLat2 = Math.cos(lat2);
        final double cosLat3 = Math.cos(lat3);
        final double lat31 = lat3 - lat1;
        final double lng31 = lng3 - lng1;
        final double lat21 = lat2 - lat1;
        final double lng21 = lng2 - lng1;
        final double a = Math.sin(lng31) * cosLat3;
        final double c = Math.sin(lng21) * cosLat2;
        final double b = Math.sin(lat31)
                + (2 * sinLat1 * cosLat3 * MathUtil.hav(lng31));
        final double d = Math.sin(lat21)
                + (2 * sinLat1 * cosLat2 * MathUtil.hav(lng21));
        final double denom = ((a * a) + (b * b)) * ((c * c) + (d * d));
        return denom <= 0 ? 1 : ((a * d) - (b * c)) / Math.sqrt(denom);
    }

    /**
     * Returns tan(latitude-at-lng3) on the great circle (lat1, lng1) to (lat2,
     * lng2). lng1==0. See http://williams.best.vwh.net/avform.htm .
     */
    private static double tanLatGC(final double lat1, final double lat2,
            final double lng2, final double lng3) {
        return ((Math.tan(lat1) * Math.sin(lng2 - lng3)) + (Math.tan(lat2) * Math
                .sin(lng3))) / Math.sin(lng2);
    }

    private PolyUtil() {
    }
}
